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1.  INTRODUCTION 


The  most  successful  discriminants  between  small  underground  explosions  and  low  magnitude 
earthquakes  are  usually  based  on  the  ratio  Pn/Lg  or  Pg/Lg  since  these  wave  trains  tend  to 
dominate  short-period  records  observed  in  continental  areas.  Pn  are  compression  waves  critically 
refracted  at  the  moho.  The  Lg  train  is  S  wave  energy  trapped  in  the  “granitic”  upper  crust.  It  can 
be  modeled  as  multiple  total  internal  reflection  of  SH  [1],  or  as  a  sum  of  higher-mode  surface 
waves  [2],  Similarly,  the  Pg  wave  train  is  trapped  P  wave  energy.  The  ratios  Pn/Lg  and  Pg/Lg 
begin  to  discriminate  between  explosions  and  earthquakes  at  frequencies  above  1Hz  and  appear 
to  work  better  with  increasing  frequency  up  to  about  9  Hz  [3], 

There  is  a  simple  rational  for  using  a  P/S  ratio  to  discriminate  between  explosions  and 
earthquakes.  An  underground  explosion  is,  to  a  first  approximation,  a  spherically  symmetric 
pressure  source  that  is  expected  to  generate  mostly  P  waves.  An  earthquake  is  a  shear  slip  along 
a  fault  plane  that  should  generate  mostly  S  waves.  The  ratio  Pn/Lg  should  be  larger  for 
explosions  than  for  earthquakes.  While  this  is  generally  true,  there  are  numerous  exceptions  that 
must  be  understood  if  this  discriminant  is  to  be  widely  used  in  different  geological  and  tectonic 
environments.  The  two  obvious  sources  of  uncertainty  in  Pn/Lg  are  source  effects  and  path 
effects.  Our  work  has  focused  on  understanding  the  source  effects.  In  specific,  we  have 
investigated  mechanisms  by  which  explosions  generate  high-frequency  P  and  S  waves  and  how 
they  depend  on  rock  type  in  the  source  region,  overburden,  and  tectonic  regime.  As  discussed  in 
more  detail  below,  our  approach  has  been  to  develop  a  micromechanical  model  for  the  dynamic 
generation  of  distributed  fracturing  (damage)  in  the  non-linear  source  region  of  underground 
explosions,  and  to  use  this  model  in  numerical  simulations  of  explosions  to  quantitatively  assess 
the  role  of  fracture  damage  in  the  generation  of  P  and  S  waves  in  the  far  field. 

Our  micromechanical  damage  mechanics  model  differs  from  the  more  commonly  used 
continuum  models  in  that  we  model  the  nucleation  growth  and  interaction  of  the  myriad  of 
fractures  driven  by  the  high  stresses  in  the  source  region.  As  we  show  in  more  detail  below,  this 
allows  us  to  physically  incorporate  variables  such  as  the  initial  fracture  distribution  in  the  source 
rock  and  loading-rate  effects  that  are  known  strongly  affect  S  wave  generation. 
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2.  BACKGROUND 


The  observed  generation  of  S  waves  by  a  spherical  explosion  source  was  first  explained  as 
resulting  from  the  dynamical  relaxation  of  the  cavity  in  a  shear  pre-stress  field  [4]  [5]  [6].  Shear 
wave  radiation  thus  generated  should  have  a  dominant  wavelength  near  the  size  of  the  explosion¬ 
generated  cavity.  More  recently,  Johnson  and  Sammis  [7]  have  shown  that  the  nucleation  and 
growth  of  fractures  in  the  source  regime  can  also  produce  significant  coherent  P  and  S  waves  in 
the  far-field.  Pre-existing  tectonic  shear  stress  and/or  a  preferential  orientation  of  the  preexisting 
fractures  that  nucleate  the  explosion-induced  fractures  enhances  the  S  energy.  The  difference 
between  this  new  “fracture-damage”  mechanism  and  the  older  relaxation  mechanism  is  that  the 
scale  of  the  fracture  sources  can  be  much  smaller  so  the  resultant  P  and  S  energy  can  be  high- 
frequency.  We  called  this  fracture-generated  energy  “secondary  radiation”  since  it  is  not 
included  in  current  seismic  source  models. 

As  detailed  in  the  following  section,  we  have  recently  improved  the  Ashby  and  Sammis  model 
[8]  by  including  high-velocity  crack  nucleation  and  growth  laws,  and  have  improved  the  Johnson 
and  Sammis  [7]  analysis  by  building  our  new  dynamic  damage  mechanics  into  the  commercial 
ABAQUS  dynamic  finite  element  code.  These  improvements  allow  us  to  properly  calculate  the 
nucleation  and  growth  of  fracture  damage  at  the  extremely  high  loading  rates  in  the  non-linear 
source  region  of  a  nuclear  explosion.  Also,  use  of  the  ABAQUS  code  allows  us  represent  the 
change  in  elastic  stiffness  associated  with  the  damage  in  a  physical  way  that  is  consistent  with 
the  damage  mechanics,  rather  than  using  the  ad  hoc  “effective  elasticity”  approach  taken  by 
Johnson  and  Sammis  [7]. 

There  is  mounting  observational  evidence  that  P/S  discriminants  are  strongly  influenced  by 
conditions  at  the  source.  Walters  et  al.  [3]  compared  the  signals  from  130  underground  nuclear 
explosions,  one  very  large  chemical  explosion,  and  50  earthquakes  ranging  in  magnitude  from  2 
to  6.  Since  the  instruments  and  propagation  paths  were  roughly  the  same  for  these  events,  they 
were  able  to  focus  on  source  effects.  For  the  explosions  they  investigated  the  effects  of 
magnitude,  depth,  and  working  point  properties  of  density,  velocity,  and  gas  porosity.  The  gas 
porosity  also  gives  some  infonnation  about  effects  associated  with  the  water  table  since  water 
saturated  rock  has  zero  gas  porosity.  Their  conclusion  was  that  Pn/Lg  and  Pg/Lg  were  more 
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sensitive  to  properties  of  the  source  medium  than  to  the  depth.  Since  velocity  and  density  are 
strongly  correlated  with  the  gas  porosity,  it  was  not  possible  to  determine  which  had  the  largest 
effect.  The  ratio  tended  to  be  larger  in  harder  rock  having  lower  gas  porosity,  with  Pg/Lg 
showing  a  larger  dependence  than  Pn/Lg.  Walters  et  al.  [3]  also  explored  the  relative  excitation 
of  Pn,  Pg,  and  Lg  at  1Hz  with  excitation  of  the  same  phase  at  6-8  Hz  for  both  explosions  and 
earthquakes.  For  all  phases,  explosions  in  media  with  high  gas  porosity  generated  relatively  less 
high  frequency  energy.  They  did  not  offer  any  physical  explanation  for  these  observations. 

Xie  and  Patton  [9]  used  a  more  fonnal  inversion  algorithm  to  separate  source  effects  from 
propagation  for  the  Pn  and  Lg  phases  from  Lop  Nor  explosions  and  nearby  earthquakes.  They 
concluded  that  “spectral  overshoot  for  explosions,  differences  in  corner  frequencies  for  Pn  and 
Lg  waves,  and  scaling  relation  for  key  spectral  parameters  are  significant  factors  controlling  the 
performance  of  the  P/Lg  discriminants  for  small  to  moderate  events  in  central  Asia.”  They  also 
comment,  “It  is  important  to  point  out  that  in  this  paper  we  made  no  attempt  to  identify  physical 
processes  such  as  cavity  rebound,  spall,  near  source  scattering,  and  interactions  between  pP  and 
P,  . . .  that  may  be  responsible  for  important  spectral  features  of  regional  phases”.  To  this  list  of 
physical  processes  we  propose  to  add  secondary  radiation  generated  during  rock  fracture  in  the 
source  and  to  investigate  its  contribution  using  our  model. 

Stevens  et  al.  [10]  considered  four  candidate  mechanisms  for  explosion  generated  Lg: 

1)  “Direct  generation  by  the  explosion  source  where  the  explosion  is  modeled  as  a  point 

compressional  source.” 

2)  “Secondary  generation  by  the  explosion  source,  where  Lg  is  generated  primarily  by  the 

nonspherical  parts  of  the  explosion  source,  with  strong  influence  from  the  free  surface.” 

3)  “Rg  scattering.  The  hypothesis  is  that  the  Rg  phase  is  scattered  as  it  travels  away  from  the 

explosion  and  is  converted  to  Lg.” 

4)  “P  scattering.  The  hypothesis  is  that  Rg  is  in  a  sense  a  variant  on  P  coda  containing 

converted  P->S  phases  in  the  crust.” 

They  favored  explanation  2.  Explanation  1  was  eliminated  because  the  spherical  source 
generates  very  little  Lg.  Explanation  3  was  eliminated  because  Rg  was  observed  to  persist  to 
large  distances.  They  found  that  inclusion  of  a  CLYD  source  with  about  half  the  strength  of  the 


3 

Approved  for  public  release;  distribution  is  unlimited. 


explosion  generates  the  observed  Lg.  This  is  approximately  the  strength  of  the  secondary 
radiation  calculated  by  Johnson  and  Sainmis  [7]  for  the  NPE  explosion.  We  propose  to  represent 
the  nonspherical  secondary  radiation  that  is  produced  by  our  model  as  an  equivalent  CLVD, 
thereby  interpreting  the  Stevens  et  al.  [10]  result  in  terms  of  physical  properties  of  the  source 
medium  that  lead  to  asymmetry  such  as  prestress,  fracture  anisotropy,  and  the  depth  dependence 
of  confining  stress,  an  effect  we  have  not  yet  explored  using  our  model. 

Our  focus  on  the  effects  of  fracture  damage  on  seismic  coupling  has  motivated  a  program  of  field 
experiments  in  a  Barre  Granite  quarry  run  by  Weston  Geophysical  and  New  England  Research 
tenned  the  Vermont/New  England  Damage  Experiment,  or  “NEDE”.  Of  particular  interest  was 
the  observed  dependence  of  both  the  induced  damage  and  the  S  wave  radiation  on  the  burn-rate 
(loading  rate)  of  the  explosive.  Also  interesting  is  the  possibility  that  the  lithostatic  gradient 
promotes  preferential  fracture  toward  the  free  surface  resulting  in  the  creation  of  a  compensated 
linear  vector  dipole  “CLVD”  contribution  to  the  source,  and  a  possible  anisotropy  in  the  dynamic 
damage  caused  by  a  known  alignment  of  the  initial  factures  (rift)  in  the  Barre  Granite  [11], 
which  is  also  a  potential  source  of  S  radiation. 

The  ultimate  goal  of  our  research  program  is  to  develop  and  verily  a  physical  model  for  the 
generation  of  fracture  damage  at  the  source  that  can  calculate  the  generation  of  secondary  high 
frequency  P  and  S  waves  for  a  variety  of  rock  types  and  tectonic  environments,  and  to  assess 
their  effect  on  current  P/S  discriminants. 

3.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

We  use  the  micromechanics  developed  by  Ashby  and  Sainmis  [8]  and  expanded  by  Deshphane 
and  Evans  [12]  which  considers  an  isotropic  elastic  solid  containing  an  array  of  penny  shaped 
cracks  all  of  radius  a  and  all  aligned  at  an  angle  ¥/to  the  largest  remote  compressive  stress  <j\ 
(Fig.l).  Sliding  on  the  array  of  penny-shaped  cracks  in  Fig.  1  produces  a  wedging  force  Fw  that 
drives  tensile  wing  cracks  to  open  in  the  direction  of  the  smallest  principal  stress  cr3  and 
propagate  parallel  to  the  largest  principal  stress  a \.  Growth  of  wing  cracks  is  enhanced  by  o\, 
retarded  by  cr3,  and  enhanced  by  a  global  interaction  that  produces  a  mean  tensile  stress  o-\" .  The 
positive  feedback  provided  by  this  tensile  interaction  stress  leads  to  a  run-away  growth  of  the 
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wing  cracks  and  ultimate  macroscopic  failure.  The  size  and  density  of  the  initial  flaws  are 
characterized  by  an  initial  damage  defined  as 


A  =  |  M  af  0) 

where  Ny  is  the  crack  density  (cracks  per  unit  volume)  and  aa  is  the  projection  of  the  crack 
radius  a  in  a  vertical  plane. 

Wing  cracks  nucleate  when  the  largest  principal  stress  reaches  a  critical  value  given  by: 


Crack  growth  is  controlled  by  sliding  on  the  cracks,  which  is  controlled  by  the  coefficient  of 
friction  //.We  ignore  the  difference  between  static  and  dynamic  friction.  We  also  assume  that 
=  45°  and,  hence,  a  =0.7.  As  the  length  of  the  wing  cracks  /  grows,  damage  increases  as 

D  =  \  Nv(l+  af  (3) 
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Figure  1:  Geometry  in  the  Micromechanical  Damage  Mechanics  Model 


Ashby  and  Sammis  [8]  show  that,  under  these  assumptions,  the  stress  intensity  factor  at  the  tips 
of  the  growing  wing  cracks  is 


i^i)(ci+c2)+  3c3 


(4) 


where 


A  - 


(5) 
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and 


c  i  =■ 


2  3/2 


(Z)/D0)1/3  1+  / 


3/2 


C2 


2  3/2 


(D!Do)'n  1 


1/2  r)2/3 


1  z> 


,2/3 


c3 


=-V 


(Z)/D0)1/3  1 


1/2 


(6) 


The  original  Ashby  and  Sainmis  formulation  [8]  is  quasistatic.  It  assumes  that  wing  cracks 
nucleate  as  soon  as  equation  (2)  is  satisfied  and  that  the  their  length  is  always  given  by  the  value 
of  /  that  makes  Kj  in  equation  (4)  equal  to  the  critical  stress  intensity  factor  Kjc  (which  is  a 
constant  for  a  given  material).  While  this  assumption  is  adequate  for  triaxial  experiments  at 
relatively  low  loading  rates  (which  were  fit  by  Ashby  and  Sainmis  [8]),  it  is  a  poor  assumption  at 
the  high  loading  rates  that  typically  occur  in  the  non-linear  source  region  of  an  explosion  (or  near 
an  earthquake  rupture  front,  or  during  an  meteorite  impact). 


When  Johnson  and  Sainmis  [7]  used  this  damage  mechanics  to  model  the  NPE  experiment,  we 
had  to  use  an  iterative  approach  in  which  we  1)  solved  the  elastic  problem,  2)  calculated  the 
equilibrium  fracture  damage  from  the  stress,  3)  adjusted  the  elastic  constants  for  the  damage,  and 
then  4)  returned  to  step  1 .  This  iteration  loop  was  repeated  until  the  damage,  stress,  and  elastic 
moduli  became  self-consistent. 


However,  at  high  loading  rates,  the  nucleation  and  growth  of  cracks  can’t  keep  pace  with  the 
rapidly  increasing  stress  field  and  true  rate  effects  become  important.  Deshpade  and  Evans  [12] 
introduced  a  simple  ad  hoc  crack  growth  law  into  the  Ashby/Sammis  [8]  damage  mechanics  to 
introduce  such  rate  effects  into  their  simulation  of  impacts  on  ceramic  armor  plating. 

3.1  Extending  Micromechanical  Damage  Mechanics  to  High  Loading  Rates 

Working  with  Professor  Ares  Rosakis  (chair  of  the  engineering  division  at  Caltech  and  an  expert 
in  high-speed  crack  growth),  we  have  replaced  Deshpande  and  Evans  [12]  ad  hoc  crack  growth 
law  with  one  based  on  the  following  theoretical  and  experimental  considerations. 
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For  a  stationary  finite  crack  under  transient  loading  conditions  the  dynamic  stress  intensity 
factor,  Ki1  evolves  with  time  following  the  application  of  loads.  Figure  2  shows  the  evolution  of 
Kjd  scaled  by  its  static  limit.  The  abscissa,  cnt/b,  denotes  the  time  from  the  beginning  of  loading 
to  the  instant  at  which  fracture  initiation  occurs  scaled  by  the  Rayleigh  wave  speed  cr  and  the 
crack  half-length  b. 

As  illustrated  in  Fig.  2,  K/1  rises  sharply  with  time,  overshoots  the  equivalent  static  value  Kst  by  a 
considerable  amount,  and  then  oscillates  around  the  static  value  with  decreasing  amplitude.  This 
oscillation  is  due  to  the  Rayleigh  waves  traveling  back  and  forth  along  the  surface  of  the  crack 
with  decreasing  intensity. 


Figure  2:  Temporal  Evolution  of  the  Dynamic  Stress  Intensity  Factor 

As  a  material  parameter,  K/c  can  only  be  obtained  through  experimental  measurements  and  is 
found  to  vary  with  loading  rate.  Figure  3a  shows  the  critical  value  of  the  stress  intensity  factor  Kj 
(scaled  by  its  steady-state  value)  for  crack  nucleation  as  a  function  of  loading  rate.  Loading  rate 
is  given  by  the  scaled  value  of  the  time  derivative  of  the  dynamic  stress  intensity  factor.  Figure 
3b  shows  the  critical  value  of  scaled  Kj  for  crack  propagation  as  a  function  of  propagation 
velocity  v. 
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Figure  3:  Normalized  Dynamic  Stress  Intensity  Factors  for  Several  Materials 

An  expression  that  represents  the  dynamic  behavior  in  Figs.  3a, b  and  Freund's  [13]  result  on  the 
dynamic  stress  intensity  factor  of  a  growing  crack  have  been  incorporated  into  a  user  developed 
material  subroutine  (VUMAT)  in  ABAQUS  to  replace  the  ad-hoc  growth  law  currently  used  by 
Deshpande  and  Evans  [12]. 


3.2  Experimental  Validation  of  Dynamic  Damage  Mechanics 

Figure  4  shows  the  failure  strength  op  as  a  function  of  loading  strain  rate  ds/dt  in  Dionysus- 
Pentelicon  marble  (from  the  Parthenon).  The  data  are  unpublished  results  of  uniaxial  and 
Hopkinson  split  bar  experiments  that  are  presented  in  our  recent  paper  [14]. 
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Figure  4:  Comparison  Between  Experimental  Measurements  and  Model  Predictions 
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This  result  shows  the  power  of  our  micromechanical  model.  Only  two  model  parameters  are 
required  (the  inital  crack  density  Do  and  flaw  size  a),  which  were  detennined  from  the 
quasistatic  uniaxial  data  at  low  strain  rates  near  10'6  s'1.  The  diamonds  are  the  model  predictions 
at  high  loading  rates,  which  clearly  capture  the  sudden  increase  in  strength  at  very  high  loading 
rates  (as  are  typical  in  the  non-linear  source  region  of  underground  nuclear  explosions). 

4.  RESULTS  AND  DISCUSSION 

Figure  5  shows  2D  simulations  of  an  explosion  using  our  new  dynamic  damage  mechanics.  A 
pressure  pulse  is  applied  to  the  surface  of  a  1  m  diameter  cavity  in  a  rock  with  initial  damage  Do 
=  0.1.  The  insets  show  the  time  history  of  this  pressure  pulse.  First  note  that  the  patterns  show 
both  circumferential  and  radial  fracturing  close  into  the  explosion  and  a  few  radial  fractures  that 
extend  to  greater  distances.  Comparison  with  fracture  maps  from  underground  nuclear 
explosions  in  Fig.  6  and  a  simulated  laboratory  explosion  in  Fig.  7  show  that  we  are,  for  the  first 
time,  correctly  modeling  the  damage  morphology.  We  simulate  explosions  in  the  lab  by 
impacting  plates  with  hypervelocity  projectiles  traveling  faster  than  the  speed  of  sound  in  the 
plate. 

Note  that  Fig.  5  also  shows  the  effect  of  the  rise  time  and  the  duration.  The  amplitude  of  the 
stress  detennines  the  extent  of  the  circumferential  cracking  while  the  duration  determines  the 
extent  of  the  radial  fractures.  Figure  8  shows  the  displacement  field  generated  by  faster  loading 
on  the  left  and  slower  loading  on  the  right.  Note  that  the  slow  loading  produces  more  radial 
fracturing  and  stronger  S  waves.  This  is  in  agreement  with  field  observations  in  the  NEDE 
quarry  experiment  which  found  that  slow  burning  explosive  produced  stronger  S  waves.  Based 
on  the  model,  the  physical  interpretation  is  that  the  S  waves  are  produced  by  asymmetry  in  the 
radial  fractures. 
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Figure  5:  Dynamic  Damage  Created  by  Different  Loading  Regimes 
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Figure  6:  Damage  Mapped  in  the  Source  of  a  Nuclear  Explosion 


Figure  7:  Damage  in  a  Glass  Plate  Produced  by  a  Hypervelocity  Impact 


12 

Approved  for  public  release;  distribution  is  unlimited. 


S  30 


§  30 


0.01 


1.00 


0.01 


1.00 


30  m 


Figure  8:  Displacement  Fields  Generated  by  Fast  Loading  (Left)  and  Slow  Loading  (Right) 

The  Ashby/Sammis  damage  mechanics  upon  which  all  this  work  is  based  is  sometimes  criticized 
because  it  assumes  that  all  the  active  flaws  have  the  same  size  and  orientation.  In  Bhat  et  al.  [15] 
we  focused  on  Westerly  Granite  (one  of  the  most  studied  rocks)  and  explored  the  effects  of 
allowing  a  range  of  flaw  sizes  and  orientations  consistent  with  petrologic  studies.  We  found  that 
relaxing  the  assumption  of  a  single  flaw  size  and  orientation  had  little  or  no  effect  on  the 
predicted  failure  surface.  The  reason  is  that  variations  in  size  and  orientation  are  compensated  by 
an  increase  in  the  density  of  flaws  that  are  required  to  fit  the  uniaxial  strength.  We  also 
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discovered  that  curvature  observed  in  the  triaxial  failure  surface  is  caused  by  plastic  flow  of  the 
weaker  minerals  in  the  granite,  and  thus  can  be  ignored  at  the  extremely  high  dynamic  stresses  in 
the  source  region  of  an  explosion. 

5.  CONCLUSIONS 

Loading-rate  plays  an  important  role  in  the  spatial  extent  and  pattern  of  fracture  damage,  and, 
consequently,  in  the  generation  of  secondary  P  and  S  waves  in  the  non-linear  source  region  of 
underground  explosions.  It  is  therefore  important  that  effects  of  loading  rate  on  the  nucleation 
and  growth  of  fracture  damage  be  included  in  simulations  of  underground  explosions.  The 
micromechanical  model  developed  here  does  this  in  a  physical  way  that  incorporates  information 
on  the  source  rock  such  as  the  size  and  density  and  anisotropy  of  the  initial  fractures  and 
information  about  the  source  environment  such  as  lithostatic  pressure  and  water  (or  ice)  in  the 
initial  cracks. 

While  we  have  built  our  new  dynamic  damage  mechanics  into  the  ABAQUS  dynamic  finite 
element  code,  it  should  be  emphasized  that  the  ultimate  goal  of  our  work  is  not  to  produce  a  new 
numerical  code  to  simulate  underground  nuclear  explosions.  Rather,  the  goal  has  been  to  validate 
the  dynamic  damage  mechanics  by  comparison  with  lab  experiments  and  identify  first-order 
consequences  for  seismic  coupling.  Our  ultimate  goal  is  to  develop  a  damage  mechanics  that 
captures  the  physics  of  S  wave  generation  and  that  can  be  incorporated  in  current  (and  more 
sophisticated)  numerical  simulators.  We  have  collaborated  with  the  relevant  private  contractors 
and  national  labs  in  the  past  to  incorporate  the  quasistatic  Ashby/Sammis  damage  mechanics  and 
plan  to  continue  these  collaborations. 
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